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Abstract 

O 

In this paper, a new rigorous numerical method to compute fundamental matrix 
solutions of non-autonomous linear differential equations with periodic coefficients is 
introduced. Decomposing the fundamental matrix solutions $(r.) by their Floquet 
normal forms, that is as product of real periodic and exponential matrices <f>(i) = 
Q(t)e Rt , one solves simultaneously for R and for the Fourier coefficients of Q via a 
fixed point argument in a suitable Banach space of rapidly decaying coefficients. As 
ry^ an application, the method is used to compute rigorously stable and unstable bundles 

of periodic orbits of vector fields. Examples are given in the context of the Lorenz 
'"" . equations and the £ 3 -model. 
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1 Introduction 

In his seminal work [1 of 1883, Gaston Floquet studied linear non-autonomous differential 
equations of the form 

V = A(t)y, (1) 

where A(t) is a r-periodic continuous matrix function of t. The main result of [T] is now 
presented, and its proof can be found for instance in [2]. 

Theorem 1.1. [Floquet, 1883] Let A(t) be a T-periodic continuous matrix function and 
denote by Q(t) a fundamental matrix solution of ([T]). Then $(£ + r) is also a fundamental 
matrix solution, $(£ + r) = $(f)$ _1 (0)$(r) 7 and there exist a real constant matrix R and 
a real nonsingular, continuously differentiable, 2t -periodic matrix function Q(t) such that 

^(t) = Q(t)e m . (2) 
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Decomposition ([2| is called a Floquet normal form for the fundamental matrix solution 
Q(t). The real time-dependent change of coordinates z = Q~ 1 (t)y transforms system into 
a linear constant coefficients system of the form z — Rz. A stability theorem demonstrates 
that the stability of the zero solution of ([!]) can be determined by the eigenvalues of the 
so-called monodromy matrix $(r). As mentioned in [2], while the stability theorem is very 
elegant, in applied problems it is usually impossible to compute the eigenvalues of the 
monodromy matrix. An even more challenging and central problem is the computation of 
the fundamental matrix solutions. The goal of the present work is to address this major 
difficulty by introducing a new rigorous numerical method to compute explicitly Floquet 
normal forms as in ^ , hence providing a direct way to obtain fundamental matrix solutions 
of ([I]). Before proceeding with a general presentation of the rigorous computational method, 
let us introduce some motivations. 

First of all, we are not aware of any method to construct rigorously Floquet normal forms 
as introduced in Theorem |1.1| Since this fundamental decomposition was introduced more 
than 125 years ago, we believe that developing a rigorous computational method leading to 
an explicit construction of Floquet normal forms is an important contribution to the field 
of differential equations. 

The second motivation is directly linked to the study of dynamical systems. Indeed, 
equations of the form arise naturally when studying stability properties of time-periodic 
solutions of differential equations y — g(y), where g : W 1 — > W 1 is a smooth map. Assume 
that r is a r-periodic orbit of y = g(y) parameterized by ^(t) S W 1 (t g [0, r]), and define 
the r-periodic matrix function A(t) — ^7g("j(t)), where V ' g is the Jacobian matrix. Consider 
$(t) the principal fundamental matrix solution of y = A(t)y — V<?(7(i))y, that is the 
unique fundamental matrix solution so that $(0) = J, an d assume that a Floquet normal 



form = Q(t)e Rt has been computed. Theorem 3.7 shows how the information from 
the Floquet normal form can directly be used to compute important dynamical properties 
of r. More explicitly, it is demonstrated that the stability of the periodic orbit T can be 
determined by the eigenvalues of R while the stable and unstable tangent bundles of T can 
be retrieved from the action of Q(t) (with t G [0, r]) on the eigenvectors of R. 

A final motivation comes from the fact that computing stable and unstable bundles of 
periodic orbits is an important step toward computing rigorous parameterization of invariant 
manifolds of periodic orbits. In fact, one of our future goal consists of combining the ideas 
of [3] to rigorously parameterize invariant manifolds of periodic orbits, and then to use 
that information to solve rigorously, following similar ideas than the ones presented in [3], 
a projected boundary value problem whose solutions would correspond to cycle-to-cycle 
connections and to point-to-cycle connections. Note that the approach of using projected 
boundary value problems to compute (non rigorously) cycle-to-cycle connections and to 
point-to-cycle connections has been adopted by several authors (e.g. see [5], [5J, [7])- 

Let us now introduce the ideas behind the rigorous method to compute Floquet nor- 
mal forms. Rather than jumping immediately into a deep mathematical description of the 
method, we present the general ideas and we refer to Section [2] for a more detailed presen- 
tation. 

The first step is to substitute the Floquet normal form $(t) = Q(t)e Rt in the differential 
equation 0. From this, it follows that (R,Q(t)) is a solution of the differential equation 
with periodic coefficients Q — A(t)Q — QR. On the converse, if a real constant matrix R 
and a 2r-periodic matrix function Q(t) solve 

Q = A(t)Q - QR . . 

0(0) = /, {i) 
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then the matrix function Q(t) :— Q(t)e Rt is the principal fundamental solution of (JlJ. 
Therefore, the problem of computing fundamental matrix solutions in the form $(<) = 
Q{t)e Rt reduces to find (R,Q(t)) satisfying The next step is to introduce a nonlinear 
operator / (see Section 2.1 for details) whose zeros are in one-to-one correspondence with the 
solutions of ([3]). Letting x = (R, Qo, Qi, Q2, ■ ■ ■ ), where the Qfe's are the Fourier coefficients 
of Q(t), the problem of computing Floquet normal forms $(t) = Q(t)e Rt is then equivalent 
to find x such that f(x) = 0. By the a priori knowledge of the smoothness of Q(i), the 
Fourier coefficients Qfc's decay fast, meaning that the solutions of f(x) = live in a suitable 
Banach space £l s of rapidly decaying coefficients. To prove existence, in a constructive way, 
of solutions of the infinite dimensional nonlinear operator equation f(x) = 0, we use rigorous 
numerics. To be more precise, the goal of rigorous numerics is to construct algorithms that 
provide an approximate solution to the problem together with precise bounds within which 
the exact solution is guaranteed to exist in the mathematically rigorous sense. It is worth 
mentioning that by now, the use of rigorous numerical methods is a standard approach to 
study differential equations and dynamical systems (e.g. see [5J, [9], [10], [11], [12], [15] . 

M- m- M- ED- 

Based on the previous discussion, the next step consists of computing a numerical ap- 
proximation x of fix) = and to demonstrate that close to x, there exists a genuine solution 
x* of f(x) = 0, corresponding to the wanted explicit Floquet normal form of the principal 
fundamental matrix solution $(i) of ([!]). However, since the operator / is infinite dimen- 
sional, a finite dimensional approximation of / must be introduced in order to compute an 
approximation x. This is done in Section [2. 2[ Once x is computed, a Newton-like operator 
T : Q, s : — > Q, s defined by T(x) = x — Af[x) is introduced, where A is an injective linear 
operator which acts as an approximation for Df(x) . Since A is injective, the fixed points 
of T and the zeros of / are in one-to-one correspondence. The next step is to consider small 
balls B s {r) C £l s centered at the numerical approximation x, and to solve for r for which 
T : B z (r) — > B s (r) is a contraction (see Section 2.3). The rigorous verification that T is a 
contraction on B s (r) is done via the use of the so-called radii polynomials which provide, in 
the context of differential equations, an efficient means of determining a domain on which 
the contraction mapping theorem is applicable. The notion of the radii polynomials was 
originally introduced in [18] and [T9] to prove existence of equilibria of PDEs. It was later 
on adapted to prove existence of equilibria of high-dimensional PDEs (e.g. see [20], [21], 
[22] , [25]h of periodic orbits of delay equations and PDEs (e.g. see [H], [25], [25], [27]) and 
connecting orbits of ODEs [I]. We refer to [25] for a more extensive and general exposure 
of the radii polynomials. 

In this work, we present a general formulation of the radii polynomials adapted to the 
context of computing rigorously Floquet normal forms (see Section 2.4 for more details). We 
present the explicit bounds in Section [2~5] that lead directly to their construction. Note that 
these bounds ensure that the truncation error terms, inevitably introduced by computing 
on a finite dimensional projection, are controlled. It is also important to mention that 
in the computation of the bounds, the floating point errors are controlled by using interval 
arithmetic .29! . In fact, all rigorous computations were performed in Matlab with the interval 
arithmetic package Intlab [30j . 

The paper is organized as follows. In Section [2] we introduce the rigorous numerical 
method to compute Floquet normal forms $(t) = Q(t)e Rt of fundamental matrix solutions 
of systems of the form ([!]) . In Section [3j we demonstrate how to use the information from 
Floquet normal forms to compute stable and unstable bundles of periodic orbits of vector 
field and how to determine the stability properties of periodic orbits. The main result of 
this section is Theorem |3.7[ Finally in Section [4] we present some applications, where we 
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construct rigorously tangent stable and unstable bundles of some periodic orbits of the 
Lorenz equations and of the £ 3 -model. 



2 Rigorous computation of Floquet normal forms 

In this section, we introduce the rigorous numerical method to compute Floquet normal 
forms = Q{t)e Rt of fundamental matrix solutions of systems of the form ([!]). As 

already mentioned in Section [T] the first step is to introduce a nonlinear operator / whose 
zeros are in one-to-one correspondence with the solutions of 

2.1 Set-up of the operator equation f(x) = 

In the following Mat(n, R), Mat(n,C) denote the space of n x n matrices respectively with 
real and complex entries. The assumption on Q(t) to be real and 2r-periodic allows to 
consider the expansion 

Q(t) = Q a + J2 Ww + ^M)^^. ( 4 ) 
fcez\{o} 

where the Fourier coefficients Qo, Qk.i € Mat(n, R) satisfy Q-k,i = Qk,i and Q-h,i — —Qk,2 
for any k > 1. Being r-pcriodic, the matrix-valued function A(t) is also 2-r-periodic, thus it 
makes sense to consider the expansion 

A(t) = ^Ae^', (5) 

fcez 

where Ao & Mat(n,R), while the matrices Ak £ Mat(n,C) satisfy A-k = C(Ak), for any 
k > 1. Here C(A) stands for the matrix whose entries are the complex conjugates of the 
entries of A. It has to be remarked that the assumption for A to be r-periodic implies that 
Ak = for k odd and A21 = Ai where Ai is the Z-th Fourier coefficient of Ait) in the basis 

{e^'jk. 

After substituting the expansions Q and ([5| in problem ([3| , the latter system of ODEs 

moves into an equation Fit) = 0, where F(t) is a 2r-periodic matrix function. By a 

k — t n, 
subsequent projection of Fit) in the Fourier basis {e 2r } ; it follows that solving (|3| is 

equivalent to solve for the unknowns 

R,Q GMat(n,R) and Q k := (Q M , Q fe , 2 ) G Mat(n,R) 2 
the infinite dimensional algebraic system 

fiR,Q ,...,Q k ,...) = (6) 

/ = fo, /l) • • • ,fk, ■ ■ ■ ) 

defined by 

/* :=Qo + 2]TQ M -/ 

fe>i 



fo := QoR - [A ■ Q)o 



(7) 



fk 



fk,l 




fk,2 





-k^Qk. 2 + Qk,iR 



' 2t 
,2tt 



{A-Q)kA 
Qk,2R - iA ■ Q) fe ,2 



k > 1 
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where (A-Q)k 1, (A-Q)k2 denote respectively the real and imaginary part of the convolution 

{A-Q) k := A kl {Q k2 ,i+iQ k2 , 2 )- 

k!+k 2 = k 

Note that /*,/ G Mat(n,R) and / fe G Mat(n,R) 2 for every fc > 1. 

The problem ([6]) consists of: ij a system of n 2 real scalar equations for /* = representing 
the initial condition Q(0) = /; ii) n 2 real scalar equations for system / = that reproduces 
< F(t), 1 >= 0; in) 2n 2 real scalar equations for each f k =Q(k> 1). Note that fk,i,fk,2 

k — t 

are the real and complex part of the equation < F(t), e 1 2t >= ik^Q k + Q k R — (A ■ Q) k . 
Here, < •, • > represents the inner product in L 2 ([0, |^]). 

Before proceeding with the analysis of the system / = given by ([6]), let us introduce 
some notation that will be adopted throughout the paper. 

Notation 

Let A, B be matrices with entries A = {aij}, B = {bi.j} and A = (Ai, . . . , A n ), 
B = (Bi, . . . , B n ) be vectors of matrices. Denote by 

i) \A\ = {\aij\} the matrix of absolute values, where | • | denotes both the real and 
complex absolute value, according with a i} j. For vectors \A\ = {\A\\, . . . , |^4,i|); 
\A\oo = ninx, ,{ a,,, ) and |^|oo = max{|Ai| oc , . . . , |^4„|oo} 

ii) A < cw B means aij < bi j for any i,j. In case b is a scalar, A < cw b means o»j < b. 
In case of vectors A < cw B and A < cw b extends as A k < C w B k and Ak < C w b, for any 
k = 1 . . .n. The same for > cw , > cw , < cw \ 

hi) | \A 1 1 oo is the standard infinity norm of a matrix: |j AW^ — max^ • \atj\; 

iv) / denotes the identity n x n matrix, 1„ is the n x n matrix whose entries are all 1. 

Coming back to the analysis of system ^ let us define the space 

x Q = (R,Q ) G Mat(n,R) 2 
x k = Qk = (Qk,i,Qkfl) G Mat{n,R) 2 ,k > 1 



X = {x = (x ,x 1 ,...,x k ,... , , ,, ,,, ny2 



Note that / : X — > X. Later on the problem of solving / = will be transformed into a fixed 
point problem for an operator T: that requires the choice of a suitable Banach subspace of 
X where to investigate the existence of solutions. To define the proper Banach space, let us 
first introduce the weigh function 

Wk -\l k = (8) 

and given x = (R, Qq, Q\ t i, Qip, ■ ■ ■ , Qk,ii Qk,2i ■ • • ) G X, let us define the s-norm of x in 
Xby 

\\x\\ s \= SUpllxfclooWfc} = SUp <^ \R\oo, |Qo|oo,SUp{|Qfc i i| 00 U'fc, \Qka\ooW s k } \. 

k>0 L k>l > 

According with the s-norm, let us define the space f2 s of sequences in X with algebraically 
decaying tails 

n s = {x G X : ||x|| s < oo}. (9) 

For any s > the space Q s endowed with the s-norm is a Banach space and the inclusion 
il s D rt s+1 holds. The introduction of n, s is motivated by the fact that a periodic solution 
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Q{t) of system ^ results to be at least as smooth as A(t). Thus, in case the function A(t) 
is analytic, it follows that Q(t) is analytic. As a consequence the Fourier coefficients of Q(t) 
decay faster than any power rate and therefore they live in Q s for any s. On the other hand, 
even a weaker assumption of the function A(t), such as a |.Afc|oo < Cj«j S for a constant C 
and positive s, allows to conclude that the solution x <E f2 s . The latter is the case we are 
mainly interested in. Indeed, for the sake of generality and to emphasize the robustness 
and versatility of the technique, one assumes the weakest assumption on Ah that makes 
the computational method applicable. Such assumption is that there exists s > 2 and a 
constant C > such that the coefficient Ak satisfy |*4fc|oo < Cw^ " . This condition implies 
that an integrable function A(t) with expansion as in ^ is differentiable up to order s — 1. 

Denote with A — {A k }k>o the sequence of the Fourier coefficients appearing in ([5| and, 
as an extension of the s-norm, define 

\\A\\ S = sup{|AU<}. (10) 

fe>0 

Lemma 2.1. Assume ||*4|| s * < oo for s* > 2. Then f maps 51 s in for any 2 < s < s* . 

Proof. Let 2 < s < s* and suppose x £ il s . Then |v4fc|oo < C\w^ s and, from Lemma 2.1 
in [20], \{A-Q) k \oo < § ■ Thus \f k {x)\ 00 <C s k\Q k \ 00 + C i \Q k \ 00 + C 2 wl s <Cw^ s+ \ for 

suitable constants C, Ci. This shows that f(x) g f2 s_1 . ■ 

Thus we will look for solutions of the system (|6| within the space W for some s > 2. 
The idea is to reformulate the zero finding problem f{x) = as a fixed point problem for a 
suitable operator T defined in VL S and to verify the hypothesis of the contraction mapping 
theorem in order to conclude about the existence of a fixed point. More explicitly, the idea 
is to prove the existence of a ball B s (r) in il s around a finite dimensional approximate 
solution x on which the operator T is a contraction. The proof will follow by verifying a 
finite number of polynomial inequalities: the so-called radii polynomials. Their computation 
will result from rigorous numerical computations and analytic estimates. The next step is 
to compute a finite dimensional approximate solution x. For this, one needs to introduce a 
finite dimensional projection of f(x) — given by 



2.2 Finite dimensional projection 

As mentioned earlier, the fist step involved in the computational method is to consider a 
finite dimensional projection and to compute an approximate numerical solution of dfjb. 

For m > 1 consider the finite dimensional space X m = JlfcLi Mat(n,M.) 2 and define the 
projections 

n m : X -> x m 

x >->n m (x) = x m = {R,Q Q ,...,Q m -\) 

IIoo : X (Qrrn Qm+l: • ■ • ) 

so that x = (x m ,H 00 (x)). Denote with 0°° := rio^O). Moreover let us define the restricted 
map 

f(m). X m X m 

x m ^n m f(x m ,Q°°) ' 

Note that for any x S X the sequence (x" 1 ,^) 00 ) € X and the finite dimensional projection 
n m applied to f{x) reads as U m f(x) — (/*, /o, . . . , f m -i)(x). Since X m is isomorphic to 



G 



K m2n , one can think of /(' 
algorithm, one computed 



an approximate zero of /^ m ' 
to identify the above vector in X 



x = (R, Q 
that is f (m \x) ? 



Suppose that using a Newton-like iterative 

■ ' Qm-l) 



0. For simplicity the same notation x is used 
and the sequence (x, 0°°) in X. As already mentioned at 
the end of Section 2.1 the idea is to consider a ball B s (r) 6 fi s centered at the approximate 
solution x and to show the existence of a contraction mapping T acting on B s (r). Hence, 
let us now introduce the fixed point operator T. 



2.3 The fixed point operator T(x) = x 

In this section, we first define an operator T on J7 S whose fixed points correspond to solu- 
tions of f{x) — and then, we introduce some computable conditions from which one can 
conclude about the existence of fixed point of T. To begin with, suppose to have chosen 
a representation of the matrices Matin, K) as vector in IR n and to have extended it to an 



to 



. Note that 



isomorphism between the space of sequences of N matrices Mat{n,] 
X rn is isomorphic to M' m2 " 2 . 

In the sequel, consider a vector V = [vi, . . . , v N2n 2 } G R N2n ■ We denote by Vk S K 2n , 
k = 0, . . . , N— 1 the vector with 2n 2 components Vk = [vk2n 2 +i, v k2n 2 +2, ■ ■ ■ , v (k+i)2n 2 }- The 
reason of this choice of notation is the following: suppose that V is the vector representation 
of the sequence x — (i?, Qq, Q\, . . . ,Qn-i) £ X N for a positive N, then Vk represents the 
couple (R,Qo) when k = and Q k = (Qk,i,Qk,2) for k > 1. 

Denote by Df^ix) the Jacobian of f^ with respect to x m evaluated at x, that is 

Df {m) ._ Df (m)^ = ^f' / ° ,/l '"""' /m ~ l) (x) e Mati2n 2 m,R). 
For clarity and completeness, 



Df 



(m) 



df* df*. 

OR OQo 

dfo dfo_ 

OR dQ 

dfx 
d(R,Qo) 



df„ 



df t df* 
dQi,i 0Qi, 2 

dfo dfo 
9Qi,i 9Qi, 2 
dfi 
dQi 



dfn 



df* 



dfo 



df, 
dQ m -i,2 

dfo 
9Q m -x, 2 



dfi 



dQ„ 



df„ 



d(R,Qo) 

where for k,j = l,...,m—l 
df k 



dQi 



d(R,Q ) 



dfk,i df k 

dR 
dfk, 2 



I dR 



dQo 
df k ,2 
dQo 



dfk 

dQ 3 



dQ n 



dfk.i df k ,i ' 



(x) (12) 



dfk, 2 



°Wj,2 
dfk, 2 



dQj,i dQj,2 . 



and each ^q'\ €E A/ai(n 2 ,IR) denotes the Jacobian matrix of the components of fk,j with 
respect to the components of Qj,i- 

Suppose to have numerically computed £)/( m ) and denote by A rn £ Mat{2n 2 rr, 
invertible numerical approximation of (D/^ 1 ™)) -1 

A m ■ DfW w I 

and for k > m, define 

A fc := &^(x) G Mat(2n 2 ,U). 



an 



dQ k 
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Lemma 2.2. Recall Q and (10 1, and assume that \\A\\ S < oo for some s>2. Then there 
exist two constants K and Ca such that for any k > K the linear operator A k is invertible 



id IIA, 



< 



The constants K and Ca depend on \\A\\ S , the period r and \R\ C 



Proof. The real and imaginary parts of (A ■ Q) k can be written explicitly as 
(A ■ Q) k ,i = {Re(Ao) + Re(A 2k ))Qks + Im(A 2k )Qk,2 + Wi 

{A ■ Q) K2 = Im{A 2k )Qk,i + (Re{Ao) - Re(A 2k ))Qk,2 + W 2 

where W\ and W 2 do not depend on Q k ,\ and Q k , 2 - Thus, looking at the definition of f k in 
Q, it follows that is of the form 



A k 



+ A 



+ A 



1,2 



2.1 



A 



2.2 



(13) 



where the entries of Ai,i and X 2t2 are linear combination of the entries of R, Ao, A 2k so that 
|Ai,i|oo, |A 2 ,2|oo < \R\oo + \Ao\oo + \A 2k \oo holds. Also, A 2 ,i = Ai, 2 only depend on A 2k . By 
a row permutation, the invertibility of A^ is equivalent to the invertibility of 



A, 



+ A 



2.1 



A 



2.2 



+ A 



1,2 



Since lAiJoo, |A 2 ,2|oo < |-R|oo + IAU + |AfeU and |A li2 |oo < lAfcU, the assumption 
\\A\\ S < oo implies that the (Ajjloo are uniformly bounded in fc, and moreover IA1.2I00 is 
decreasing. Thus there exists K such that A^ is diagonally dominant for any k > K. This 
is enough to conclude that A^, is invertible for any k > K. 

Denote by atj the diagonal elements of A k . Hence, if A k is diagonally dominant, that is 



if Kil > J2j^i \A k (i,j)\ for any 
following bound 



1. 



< max 



, 2n 2 , then using a result from 
1 



], one gets the 



Therefore, for k > K 

||oC = I. M x 

for a constant Ca depending on r, R, |^4q|oo and |^2fc|oo- 




Suppose that we chose the finite dimensional parameter to > K where K, as defined in 
Lemma [2.2| is such that A^ is invertible for any k > K. A formal diagonal concatenation 
of the operator A m and the sequence ATT , for k > to, produces the linear operator 



{Ax) k := | 

We define the operator T on X as 



A : X X 
x i— > Ax 

A m x m ) k = 0, . 

-l k 
A k Q k k > to. 



(14) 



T(x) 



Af(x), 



and denote T k (x) = (T(x)) k 
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Lemma 2.3. Recall Q and (10), and assume that \\A\\ S * < oo for s* > 2. Then for any 
2 < s < s* , T : n s — > fl s and solutions ofT{x) — x correspond to solutions of f(x) = 0. 



2.1 



given x £ fl s it follows that f(x) £ tt s 1 . The linear operator A 



Proof. From Lemma 

maps f2 s_1 in Q s . Indeed for k > m, \(Ax)k\oo = |A / 7 1 Xfe 

and assuming x £ f2 s_1 , it follows that |(Ac)fc| 00 < 

j, CI. This proves that T : ft s ' — > £l s . Since A m is invertible by assumption 



Lemma 



2.2 



< Ha; 

c\\i 



•&k | oo ■ 

Ci 



< 



Thus from 
for positive 



constants 

and Afe have been proved in Lemma [2. 2| to be invertible for all k > m > K, it follows that 
the linear operator A is invertible and therefore fixed points of T correspond to zeros of 
fix). ■ 

By construction, when restricted to the finite dimensional reduction II m 57 s , the operator 
T acts as T(x m ) = x m - A m f {m \x m ). Thus on II m ft s , T is close to the Newton operator: 
the only difference is that point where the derivative is computed does not change along 
the iteration process. Therefore we can consider T as an extension to a infinite dimensional 
space of a finite dimensional Newton-like operator. 

The existence of a fixed point for the operator T will be assured by the Banach Fixed 
Point Theorem once the operator T has been proved to be a contraction on a suitable ball 
in f2 s . The suitable ball on which T will be proved to be a contraction will be sought within 
the family of balls B s (r) £ Q, s 

B s (r)=x + B(r) 

where B{r) is the ball of radius r in fi s centered in the origin and r is treated as variable. 
Following the same approach as in different other papers (e.g. see [20], [16], [24], [26] . 
[27] . [25] . [22], [23], [4], [18], [H]), we are going to construct a finite set of computable 
conditions, the so-called radii polynomials, to be solved in r, whose verification implies that 
the hypothesis of the Banach Fixed Point Theorem are satisfied. In practice, the radii 
polynomials are defined as realization of the hypotheses of the following theorem. 
Suppose there exist two matrices sequences 



Y = (Y , Y u ... Y k) . . . ), Z(r) = (Z Q , Z u ... Z k , . . . )(r), 



Y,Z£X 



such that 



\(T(x)-x) k \ < c 



sup 

6i,6 2 eS(r) 



[DT(x + 6i)6 2 ] k < cw Z k {r), Vfc > 0. 



(15) 



Theorem 2.4. Fix s > 2 and let Y and Z defined as in (15). If there exists r > such 



that || Y + Z \\ s < r, then the operator T maps B s (r) into itself and T : B s {r) — > B s (r) is a 
contraction. Thus, by the Banach Fixed Point Theorem, there exists an unique x* £ B s (r) 
solution ofT(x*) — x* and therefore solution of f{x*) = 0. 

Proof. Two statements need to be proved: 

i) T(B s (r)) C B s (r), that is \\T(x) - x\\ s < r for all x £ B s (r), 

ii) T is a contraction, that is there exists k £ (0, 1) such that for every x,y £ B s (r), one 
has that \\T(x) - T{y)\\ s < k\\x - y\\ s . 

For a given k > and any x,y £ B s {r), the mean value theorem implies 

T k (x)-T k (y) = DT k (z){x-y) 
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for some z £ {tx + (1 - t)y : t G [0, 1]} C B s (r). Note that r ife^r € B(r) thus for ([15]) 

r(x - y) 



\T k (x)-T k (y)\ = 



DT k (z), 



f - y\\ 



-\\x - y\\ s < cw Zk T \x - y\ s (16) 
r r 



The triangular inequality applied component-wise gives 

\T k (x) ^ %k\ <cw \T k (x) - T k (x)\ + \T k (x) - x k \ < cw Y k + Z k (r) 

hence 



\T k (x) - Xfcloo < \Y k + Z k (r)\ c 



Therefore for any x G B s {r) 



\\T(x) - ill, = sup{|T fe (x) - < sup{|r fe + Z fc (r)U<} = \\Y + Z(r)\\ s < r. 

k>0 k>0 

This proves i). 

Again from pi, for any x,y e B s (r), \T k {x) - T^y)^ < ^ Zk ^°° \\x - y\\ a , thus 



\\T(x)-T(y)\\ s <^P^\\x-y\\ 



(17) 



Note that all the entries of Y k and Z k (r) are non negative, thus 1 2^ (r) | oo < \Y k + Z k (r) 
and \\Z(r)\\ s < \\Y + Z(r)\\ s < r. That implies that 



«:=M* e(0il) , 
r 

and we can conclude the proof of ii). An application of the Banach Fixed Point Theorem 
on the Banach space B s (r) gives the existence and unicity of a solution x* of the equation 
T(x) — x in B s (r) and, from Lemma 2.3 of a solution of f(x) =0. ■ 

2.4 The radii polynomials 

As already mentioned in Section [T] the radii polynomials are a set of r-dependent polyno- 
mials p k {r) defined in such a way that if r* is a common solution of p k {r*) < 0, then the 
ball B s (r* ) C f2 s of radius r* centered at the numerical approximation r* contains a unique 
solution of f(x) — 0. This is due to the fact that by construction of the polynomials, one 
has that ||Y + Z(r*)\\ s < r* , meaning that the hypotheses of Theorem 2.4 are satisfied. In 
terms of the components, the formula || Y + Z(r) || s < r reads as 



\Y k + Z k {r)\ c 



< 0, Vfe > 0. 



w. 



(18) 



The latter consists of a system of infinitely many inequalities, which is then impossible to be 
verify directly with computations. In order to reduce ( 18 1 to a finite number of inequalities, 
suppose that, for a given M, there exist Ym and Zj^[r) such that 



\{T{x) - x) fe | 00 < 



M s 



-Y, 



M - 



sup 

b 1 ,b 2 £B(r) 



[DT{x + b 1 )b 2 ] 



M s 



< ~Z M (r), Vfc > M, (19) 

oo rC 



and introduce the set of M + 1 radii polynomials as follows. 
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Definition 2.5. The radii polynomials are defined as 

Pk(r) := Y k + Z k {r) - ^(1„, 1 B ), k = 0, . 
Pm ■= Y M + Z M - 7zs-- 



,M 



(20) 



Theorem 2.6. Consider M and let Y, Z such that Y k , Z k satisfy ( 15 ) for k = 0, . . . , M — 1 

while for k > M define 

M s M s 
Y k := — r M [l„,l„], Z k (r) := — Z M [t n , l n ], 

where Ym,Zm satisfy the tail condition (19). If there exists r > suc/i that p k (r) < cw /or 
all k = . . . , M, then there exists a unique x* £ B s (r) such that T(x*) — x* and f(x*) = 0. 

Proof. Since by definition Y k > cw 0, Z k > cw 0, the relations Pfc(r) < clu imply that |Yfe + 
Z k (r)\oo < for jfe = 0, . . . ,M- 1. For fc > M, F fe , Z fc satisfy and from Yjvf + Zm. (r) - 
^ < 0, it follows that |Y fe + ^(rJloo - ^ < 0. Indeed 



Hence 



\Y k + Z k (r)\ 00 = *f(Y M +Z M )<^^ Vfc>M. 



||F + 2|| s = sup{|F fe + ^K} <r, 

fe>0 



and the result follows from Theorem 12.4 



2.5 Construction of the bounds Y, Z 



This section is devoted to the construction of the matrices Y k , Z k satisfying (15), and of 
the asymptotic bounds Ym, Zm satisfying (19). This construction provides the complete 



description of the radii polynomials introduced in Definition 2.5 With the aim of remaining 
as general as possible, the only constraint we assume on the r-periodic function A(t) is that 
the vector of Fourier coefficients A given in ^ satisfies ||^4|| s * < oo for s* > 2. Nevertheless, 
further information on the coefficients A k may be useful to get sharper analytical estimates. 

In what follows, the growth rate parameter s has been fixed so that 2 < s < s*, the 
finite dimensional parameter m has been chosen greater than K, where K is a lower bound 
given by Lemma [272] and the computational parameter M has been chosen so that M > m. 
Moreover, assume that one computed A^T 1 for k = m, . . . , M — 1. Note that in some cases, 
it will be possible to achieve this task analytically, but in other cases, only an interval 
enclosure using rigorous numerics will be possible. Also, recalling Lemma [2~2| denote by Ca 
a computable constant such that 



< 



C A 



for k > 



(21) 



2.5.1 The bound Y 

By definition, T(x) — x — —Af(x), thus define Y as 



Y k 



{A m fW(x)) k \ 



k = . . . , m — 1 
k = m, . . . ,M - 1. 



(22) 
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The tail bound Y M is defined so that Ym^t > | A^ 1 f k (x) | , for any k > M. We now 
introduce a coarse bound Ym based on the relation |A^ 1 /fe(x)| < ||A i T 1 || co |/ / ! C (x)| 00 . Since 
Qk,i — Qk,2 — for any k > to, it follows that 



fk{x) = 



-(A-Q) ktl 

~(A ■ Q)k,2 



E 

k%+k 2 =k 
\k 2 \<m 



-Re (A kl {Qk 2 ,i + iQk 2 a)) 
-Im (A kl (Qk 2 .i + iQk 2 a)) 



Vfc > M. (23) 



Now, using the fact that |-4fc|oo < ||.4|| s *io fe a , both \f k: i(x)\ and \f k ^{x)\ are component- 
wise bounded by 



Ak 1 (Qk 2 ,l + iQk 2 ,2) 



ki+k 2 =k 
\k 2 \<m 



— CW 

^WQk,,! + iQk 2 a\ 

k%+k 2 =k 
\k 2 \<m 

m— 1 

<cw \A k \\Qo\ + d-^fc-'l + l^fc+ll)IOi,l + *Ql,2\ 



1=1 



m—1 



-i„|Qo| + $>j 



i=i 



u k-i 



K\Qi,i +iQi,2\ 



For k > M the bounds 



< 1 and w 



i 



< 1 + (l — jj) hold, thus one 



computes the matrix 

W = t n \Q \ + J2 1^+ I 1 

so that 

|A(aO|oo < fc" s 
Finally, using || A^" 1 < define 



m—1 



M 



J t n \Qi,i + iQia\ 

.|W|oo, forfc>M. 



Y M := 



-MIU^aIW^U. 



2.5.2 The bound Z 

To construct the bound Z so that 

sup [DT(x + b 1 )b 2 ] k < cw Z k (r), Vfc>0, 

bi,b 2 £B(r) 

it is convenient to factor the points 61,62 € -B(r) as 61 = ™, 62 = rw with u, v £ -B(l), to 
expand in the variable r and finally to uniformly bound the expression using the fact that 
u, v e B(l). Denote u = [u , Mi, ... , w^, . . . ], where each u k — (ufe,i, Ufe,2) € Mat(n, M) 2 . In 
order to simplify the exposition, both the matrices Wfc,i,Ufc,2 will be denoted as u k . Indeed, 
what really matters is the bound \u k ,2\ < C w w k s that finally will be applied to obtain 

the uniform estimates. The similar notation for v k . 

Let us introduce the linear operator : Q s+1 — > s defined as 



{A t x) [ {Df^. X ™) k k = 0,...,m-l 
( AfeXfe, k>m, 



(24) 
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and consider the splitting 



DT(x + ru)rv = [I — ADf(x + ru)] rv 

= [I - AA^] rv- A [Df(x + ru) - A>] rv. 



(25) 



The definition of Z(r) will follow as a result of different intermediate estimates: indeed we 
are going to introduce the vectors Z°, Z 1 , Z 2 such that 

\[I - AA^]rv\< cw Z°r , Vv G B(l), 
I [Df(x + ru) — A^]rv\ < cw Z x r + ZV, Vu,«€-B(l). 



From (25) it follows that 



[DT(x + ru)rv] 



[(I - AA*)rv] k + [A(Df(x + ru)-A^)rv] k 



Hence, the elements Z k , for k = 0, . . . M — 1 can be defined as 

|4 ffl |(^r + Z 2 r 2 ) m 
7j l 



Z k (r) = > Z * r 



Z° k r 



K\Zir + Zlr 2 ), 



k = 0, . . . , to — 1 
k = to M — 1. 



Finally the element Zm will be defined to satisfy (19). 



The bound Zq 
Since \v k \ < cw W /T s l> define Z° as 



(z o } f [\I-AAt\{wJ'i%£] 



k = 0,. 

k > m 



, TO — 1 



(26) 



(27) 



so that 

\[I-AA^\ rv\ < cw Z°r. 

Note that A 1 is an almost inverse of A, indeed by definition A m D/( m ) « /. Then the size 
of Z° is small and depends on the accuracy of the numerical method that computes the 
inverse A m . 

The bounds Z 1 , Z 2 

Concerning the terms in Z 1 , Z 2 , consider the expansion as quadratic polynomial in r 



[(Df(x + ru)-A^)rv] k = £ c M 



(28) 



First note that 



(A-Q) 



fc,i 



(^Re(Aj)Qi.i - Im{Aj)Qi,2 



j+i=k 



(A-Q) k ,2= (l m (Aj)Qi,i + Re(Aj)Qi,2 i 
j+i=k 

then, taking in mind that Qfe.2 = —Q-k,2 and denoting with sg(l) = sign(Z), one computes 



co,i = 



(Re{A 3 ) - sg{l)Im{A 3 ))v w 



l+j=0 
\l\>m 



C0 : 2 





UqVq + VqUq 



(29) 
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for k = 1, . . . , m — 1 



Ck.l 



E 

i+ 3 =k 

\l\>m 



and for fc > m 



Ck.l 



E 



fi?e(^) - s ff (/)Jm(^j))"|i| 
(lm{Aj) + sg{l)Re(Aj))v m 



Re(Aj) - S5f(/)7m(^))u|;| 



Cfc,2 



(30) 



Im(Aj) + sg(l)Re(Aj) 



"W 



Cfc,2 



UkV + V k U 

UkVo + v k u 



(31) 



Therefore Z 1 , Z 2 need to be defined so that Z\ > cw \ck,i\ and Z 2 > cw \ckfi\- To achieve 
this, it is enough to substitute in the above expression the bounds |ufe|, \vk\ < C w w k ~ s t n and 
| ± Re(Aj) ± Im(Aj)\ < cw \Re(Aj)\ + \Im(Aj)\. Since 1„1„ = nl„, one gets 



|co,2| <cw 2n 



\ckM < cw 2nw k ' 



71 



=: Z fc 2 , k > 1. 



(32) 



However this approach is not completely feasible for the computation of \ck,i\, due to the 
presence of series. Therefore it is necessary to introduce further computational parameters 



L k > max{fc, to} 



(33) 



and matrices H k so that 



I c , 1 1 <c 



|Cfc.l| <c 



2 (|J?e(^)| + |7m(^)|)^- s lr 



m<|/|<L 



Ho =: Z n 



E 



m<\l\^k,\l\<L k 



(\Re(A 3 )\ + \Im(A 3 )\)w^t n 
(|7TO(^)| + |i?e(^)|) W pl„ 



, TO — 1 



(34) 

and similarly for k > m. It means that the bound Z 1 has been defined as sum of two 
factors: the first obtained by rigorous computation of a finite number of elements in the 
series, the second analytically defined to estimate the tail part of the series that have not 
been computed. 
Define 



C(M,s) :-- 



1 



+ 



1 



+ 



(M + l) s (M + 2) s s - 1 (M + 2)*- 1 ' 



and 



where for k > 



H n := 



2((L ,s)t n 
h t n 



H k :— h k 



1„ 
1„ 



(35) 



hk 



V2n\\A\l 



(L k + 1- ky 
Hence, one has the following result. 



-(((L k -k,2 S ) + ((L k ,2s) 
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Lemma 2.7. Formula (34) holds for H , defined in (35 1. 



Proof. First note that for any M > 1 and s > 2 

oo 

E ri<C(M, S ). 



fc=M+l 



(36) 



That can be seen from the fact that X)a*La/+i = (m+i) s 



(M+l) = ^ (M+2Y 



Eoo 
fc= 



fc=M+3 fe s 



< 



(M^y + (m+2)° + Im+2 ^ Hence, one has that 

oo 

For the remaining terms note that (\Re(Aj)\ + \Im(Aj\) < cw V%\Aj\ < cw V2 ',i* l ra , then, 

3 

for any k > 0, the tail part of the series can be bounded by 



]r (mA^ + \i m (A 3 )\)wr t n 

V2n\\A\\ s * 



i+j=k 

\l\>L k +l 



°° ( 1 1 \ 



(L fc + 1 - ky 



v r 1 i \ W - S1 

^—^ \ in? . in? . / 1 



l=L k + l V l ~ k k + l ■ 

In the last passage we have used the fact that > k, s* > s and the relation l rl l„ = nt n . 



The result follows by applying ( 36 ) once the last series has be rewritten as 

y (— + — 

— ^ \ in? . in? 



l=L k +l 



■ J i-k uy fe^ 



E ( ,„? 



l=L k + l ^ W k+l 

OO 

< E -r 2 

l=L k +l 



E 



l=L k -k+l V 



E 



w, 



l=L k -k+l 



The bound Zm 
From [23], one has that 



V — -< — 2 + 2) - 



- h w k™k~ w 



k 1 +k 2 =k K i K 2 
|fci|#* 



where 









Vk = 2 




+ 



41og(Jfe-2) tt 2 -6 



+ i?m - 1 - 



w. 



2 k 



2 1 

k + 2 
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Recall the definition of cj^i and c k ,2 given in (28 1, with a more explicit form in (31 ) for the 
case k > m. Then for k > M one has that 



Cfel 



<V2\\a\\, 

V2\\A\\ 



• E 

l+j=k 
\l\?k 



— s —s 

w j w i 



< 



Cfe,2 



< 



2n 



M 



1=1 



Vm - 1 



Since for fc > M the first term on the right hand side of (26) is zero, the following estimate 
holds 



[DT(x + ru)rv] k 



< IIA 



[A (Df(x + ru) - A*)rv] k 
fe 1 ||oc(|cfe,i| 00 r + \c kt2 \<x,r 2 ). 



(37) 



Finally, combining (21) and that k > M, one gets that ||A fe ||oo 



< and we can define 



Zm as 



Z 



M 



Ca 
M 2 



(V2\\A\\ s *C 1 r + 2nr 2 ). 



3 Computing stable and unstable tangent bundles of 
periodic orbits using Floquet normal forms 



Consider an autonomous differential equation 

y = g(y), geC^R") 



(38) 



and suppose that 7(f) is a r-periodic solution with 7(0) = 70. Denote by T = {7(f), f £ 
[0,t]} the support of 7 and for any 6 e [0,r], define 79 (f) = 7(f + 9) the phase-shift re- 
parametrization of T. Being autonomous, system (38) has the property that any of the 
curves Jg(t) is a r-periodic solution satisfying jg (0) = 7(0) . We refer to T as the periodic 
orbit and 79 as the periodic solutions. 



Definition 3.1 (Monodromy matrix). Let 7 : R — >• K™ be a r-periodic solution (38) and 
let <&e(t) be the unique solution of the non-autonomous linear problem 



*e = Vg(7e(i))$ e 
$9(0) = /. 



The matrix $e(r) is called the monodromy matrix of 70(f). 



(39) 



Having chosen 7(f) = 70(f), in the following we identify $(t) = $o(' r )- The next two 
Lemmas are classical results and are direct consequence of <$>g(t) being a fundamental matrix 
solution. For sake of completeness, we present their proofs. 



Lemma 3.2. For any 6 € [0,r], f/ie solution $e(f) 0/ (38) satisfies 
<t> g (nT + f) = $ e (f)$e(r) n , Vf eR, Vn£ N. 
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Proof. Without loss of generality let us consider 6 = 0. By induction on n > 0. For n = 
the result is obvious. Suppose it holds for n — 1. Then 

$(nr) = $((n - l)r + r) = ^(t^t)™- 1 = $(r) n . 

Define 

= $(* + 7M-)$(tm-) -1 . 

It follows that *(0) = J and that 

= $(nr + t)$(nr) _1 = A(nr + £)<J>(nr + i)$(nr) _1 = A(f)*(t) 
For the uniqueness of solutions of the initial value problem, ^f(t) = thus 
$(f + nr) = $(t)$(nr) = $(t)$(r) n , Vf G E. 

■ 

Lemma 3.3. The matrices ^g(r) are equivalent under conjugation. In particular 

$ 9 (r) = $(0)$(t)$(0) _1 . (40) 

Proof. The matrix $(t) := &(t+6) is solution of the equation y = Vg(j(t+0))y = Vg(7e(t)), 
with $(0) = $(0). Since is the principal fundamental solution of the previous system, 

*(t) = 

It follows 

*»(t) = $(i)$(6')- 1 = *(t + 0)$(0)-\ V*. (41) 

Thus 

$ (r) = $(r + 0)$(0)~ 1 = $(6»)$(r)$(6»)- 1 
where, in the last passage, Lemma [3. 2| has been used. ■ 



The previous result implies that all monodromy matrices <S>g(r) have the same eigenval- 
ues. That motivates the following definition. 

Definition 3.4. The eigenvalues Oj of the monodromy matrix $(r) are called the Floquet 
multipliers of the periodic orbit T. 

As already mentioned in Section [T] in the theory of dynamical systems the monodromy 
matrix $(t) associated to a periodic solution j(t) plays a fundamental role since it encom- 
passes the information about the stability character of 7. Indeed, as shown in Proposi- 
tion 2.122 in [2], the Floquet multipliers of j(t) are in fact the eigenvalues of DP(j(0)), 
where P(x) denotes the Poincare map of j(t) on a (n — 1) -dimensional hypersurface transver- 
sal to 7 at 7(0). Moreover, it can be proved that at least one of the Floquet multipliers aj 
of $(t) is equal to one, corresponding to the eigenvector 7(0). Hence, we will denote by 
cr n = 1 the Floquet multiplier corresponding to 7(0) and denote by {o-j}j=i ... n -i the set of 
non trivial Floquet multipliers. We refer to Section 2.4 in [5] for a more extensive analysis 
of the links between Poincare sections and Floquet theory. Based on the above discussion, 
we are now ready to introduce the definition of stability of a periodic orbit. 



Definition 3.5. Let V = {~/{t),t <G [0, r]} be a r-periodic orbit of the system (38) and let 
{oj}j=i n— 1 be the corresponding set of non trivial Floquet multipliers. We say that 
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• r is stable if V j € {1, . . . , n — 1}, \oj \ < 1; 

• r is unstable if 3 j € {1, . . . ,n — 1} such that |<Tj| > 1. 

Moreover, if p < n — 1 Floquet multipliers have modulus less than one, and q < n — p 
Floquet multipliers have modulus greater than one, T is said to have p stable directions and 
q unstable directions. 

Let us mention that there is a variant to the Floquet normal form introduce in The- 
orem namely there exist a constant (possibly complex) matrix B and a nonsingular 
(possibly complex) continuously differentiable, r-periodic matrix function P{t) such that 
$(*) = P{t)e Bt . We refer to Theorem 2.83 in [2] for more details and for the proof. There- 
fore, there exists a (possibly complex) matrix B such that $(r) = e Br . Denoting by Xj the 
eigenvalues of B, it follows that aj — e TXj is a Floquet multiplier. Note that for a given aj, 
the solution Aj of aj = e rXj is not uniquely defined. Indeed for any k G Z, e T ^ Xi+lM ^ L ^ 1 = aj. 
This reflects the fact that in the complex Floquet normal form <f>(i) = P(t)e Bt , the matrix 
B is also not uniquely defined. In the literature it is common to call a Floquet exponent 
associated to aj any complex number Aj so that aj = e ■» . On the converse, for any aj 
there is a unique real number lj so that \aj \ — e ljT . That motivates the following definition. 

Definition 3.6. A Lyapunov exponent associated to a Floquet multiplier Oj is the unique 
real number lj so that \aj\ = e ljT . 

Note that using the notion Lyapunov exponents, a definition of stability of a periodic 
orbit similar to the one of Definition |3 . 5| can be introduced. Indeed, given a r-periodic orbit 



r = {j(t),t € [0, r]} of (38) and considering {lj}j=i,...,n-i to be the corresponding set of 
non trivial Lyapunov exponents, we say that L is stable if L < 0, V j — 1, ... ,n — 1 and 
that r is unstable if there exists j G {1, . . . ,n — 1} such that lj > 0. 

Given a real n x n diagonalizable matrix A, let us introduce the notation = 
{ctk,Vk}k=i,...,n to denote the eigendecomposition of the square matrix A, i.e. Av^ = afcVfe, 
for all k = 1, . . . ,n. 

The following result shows how the information from the couple (R,Q(t)) coming from 
the Floquet normal form $(t) = Q(t)e Rt can directly be used to study the dynamical 
properties of the periodic orbit T. More explicitly, it demonstrates that the stability of T 
can be determined by the eigenvalues of R while the stable and unstable tangent bundles of 
r can be retrieved from the action of Q{t) (with t € [0, r]) on the eigenvectors of R. 



Theorem 3.7. Assume that T — {^(t), t £ [0,r]} is a T-periodic orbit of (38) and consider 
&(t) the fundamental matrix solution of the non- autonomous linear equation y = Vg(^(t))y 
such that $(0) = /. Suppose that a Floquet normal form decomposition of Theorem 
$>(t) = Q(t)e Rt is known. Assume that the real n x n matrix R is diagonalizable and let 
E(i?) = {fJ-j, Vj}j=i,...,n the eigendecomposition of R. Then the Lyapunov exponents lj ofT 
are given by 

lj = Re{tx 3 ). (42) 
Furthermore, for any 6 G [0, r], if one defines 

w°j :=Q(0)vj, (43) 

then uij is an eigenvector of &e( T ) associated to the Lyapunov exponent lj. Note that Wj is 
a smooth 2r-periodic function of 6. 
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Proof. Consider the eigendecomposition = {/J-j, fj}j=i n °f the diagonalizable ma- 

trix i?, meaning that the set {vi, . . . , v n } consists of n linearly independent eigenvectors of 



R. By Lemma |3.2| one has that $(t) 2 = $(2t). Since Q(t) is 2-r-periodic and Q(0) = /, 
it follows that <1>(2t) — e R2r . Since R is diagonalizable, $(2r) = e^ 2 " 1 " is also diagonaliz- 
able. Since $(2r) = $(t) 2 and since the matrix $(r) is invertible and defined over the field 
of complex number (which has zero characteristic), then it can then be showed that 5>(t) 
is also diagonalizable. Now, since $(2r) = $(t) 2 one has that if (<r, w) £ S($(r)), then 
(er 2 , w) £ E($(2r)). Combining this last point with $(r), $(2r) being diagonalizable implies 
that the eigenspaces of $(t) and $(2r) are in one-to-one correspondence. That implies the 
existence of a set {Gj}j=i,,..,n such that S($(r)) = {crj, Vj}j = i t ,„ l7l . From the property of the 
exponential matrix operator, E($(2r)) = {e Mj2r , i>j}j=i,... in = E($(r) 2 ) = {cr 2 , Wj}j=i,...,n. 
This implies that er 2 = e Mj2r for any j = 1, . . . , n. Note that ij = Re(nj) is the unique real 
number so that |<r 3 '| = e' iT . Hence, lj is a Lyapunov exponent associated to the Floquet 
multipliers (jj. 



Now, from pT[ ), one has that 

$ 9 (2r) = $(2t + 0)$(0) _1 = Q^e^V^Qfff)" 1 , We [0,t] 



thus 



<f>e(lT)Q(9)vj = Q(9)e 2r \ = e 2T ^Q(9) Vj 



showing that E($g(2r)) = {e 2rAlj , Q(6')wj}. Applying the same argument than above, one 
can conclude that £(<J>£i(t)) = {oj, Q(#)wj}j=i,..., n forms an eigendecomposition of the ma- 
trix ^g(r). Hence, u>| = Q(8)vj is an eigenvector of $e(r). By the smoothness and the 
2r-periodicity of the matrix function Q(9), one can conclude that u>| = Q(9)vj is also a 
smooth 2r-periodic function of 9. ■ 



Recall (43) and consider w® = a® + ibj. We define the stable and unstable subspaces 



E° s , El C T 7(e) E™ of the periodic orbit T at the point 7(6*) as 

^=5j»n{o?,6?:|(7 i |<0} 

252 = SpanK, if :| ffi |>0}. 
That allows us to define the following 

Definition 3.8. We define the stable and unstable tangent bundles of T respectively by 

E S ,E U c T r R n 

E s = |J {7W}x£ s e , E u = |J {7((?)}x< 
ee[o,T] 9e[o,r] 



It is important to remark that from the conclusion of Theorem |3.7[ the complete structure 
of the stable and unstable bundles can be recovered by the action of the matrix function 
Q{i) on the eigenvectors of R, which themselves correspond to the stable and unstable 



directions at the point 7(0) on T. Also, the proof of Theorem 3.7 is constructive in the 
sense that combined with the rigorous computational method of Section [2] it provides a 
computationally efficient direct way to obtain the eigenvectors W® of $g(r), which are the 
ingredients defining the bundles of Definition |3.8| Note that one could be tempted to 
use the fact that $(r) = Q(r)e- Rr and then attempt to compute the eigendecomposition 
of $(t) directly. However, that would imply having to compute the exponential of an 
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interval valued matrix, which turns out to be a difficult task (e.g. see [32], [33]). This 
being said, the rigorous computation of the eigendecomposition of the interval matrix R is 
not completely straightforward. We addressed this problem by adapting the computational 
method based on the radii polynomials in order to enclose all the solution {/life, Vfc} of the 
nonlinear problem (R— nl)v = with constrain |w| 2 = 1. Further details on the enclosure 
of the eigendecomposition of interval matrices are postponed on a future work of the same 
authors [31]. 



4 Applications 

In this section, we present some applications, where we construct rigorously tangent stable 



and unstable bundles of some periodic orbits of the Lorenz equations in Section 4.1 and of 
the £ 3 -model in Section 4.2 Note that all rigorous computations were performed in Matlab 



with the interval arithmetic package Intlab |30j . 

4.1 Bundles of periodic orbits in the Lorenz equations 

Consider the following three dimensional system of ODEs, known as the Lorenz equations 

111 = c(u 2 - ui) 

11 2 = pui — u 2 — uiu 3 (44) 
u 3 = uiu 2 - f3u 3 

with the classical choice of parameters f3 = 8/3, a = 10 and p left as a bifurcation pa- 
rameter. Suppose to have rigorously proved the existence of a real r 7 -periodic solution 



l(t) = [7 ,7 ,7 ](*) of ([44]) in the for 

ke 



7 J W=E^ fc ^' J = 1-2,3 (45) 



in a ball of radius r 7 and centered at [r 7 ,£fe], |fc| < M 7 , with respect to the Q s norm, 
meaning that 

\Ty-fj\ < T 1 

\Re{i k )-Re{l k )\^<r^wl s \ |Jm(&) - /m(£ fe )U < r^ s \ |fc| = 0, . . . , M 7 (46) 

l-Refe)^ < rvV**, llm^k)]^ <r jW - s \ \k\ > M 7 

for a decay rate s* > 2. Note that ^ € Z 3 and ^_fe = £(£&). The existence of such 
solution could be achieved by applying a modified version of the method discussed in the 
previous section. Even with some technical differences, the philosophy is the same. Rewrite 
the system of ODEs as a infinite dimensional algebraic system where r 7 and the Fourier 
coefficients are the unknowns, then consider a finite dimensional projection and compute 
a numerical approximate solution f 7 , {£k)k- Then, by means of the radii polynomials, prove 
the existence, in a suitable Banach space, of a genuine solution r 7 , (£&)& of the infinite 
dimensional problem in a small ball containing the approximate solution. Note that this is 
not the first time that the radii polynomials are used to prove existence of periodic solutions 
of differential equations (e.g. see Q£], [H], [2B], [27], [25]). 

In the following we aim to combine the rigorous computational method of Section [2] 



together with Theorem 3.7 to rigorously compute the stable and unstable tangent bundles of 
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the periodic orbit j(t) given by (45). This first requires the computation of the fundamental 
matrix solution of the linearized system along j(t), that is the solution for t £ [0, r 7 ] of the 
non-autonomous system 

' $ = V$(7(*))* 
$(0) = I 



(47) 



where g is the right hand side of (44), Vg denotes the Jacobian of the right hand side of 
system (44) and / is the 3x3 identity matrix. The former system is nothing more than a 
particular case of ([!]), where A(t) = Vg(7(t)) and n = 3. We now apply the computational 
method presented in Section [2] to compute the principal fundamental matrix solution of the 
non-autonomous linear system y = V g(j(t))y. In particular a constant matrix R and the 
Fourier coefficients Qk of a 2r 7 -periodic function Q(t) will be computed, so that 

$(*) = Q(t)e Rt 

is the unique solution of ([47] ) . Once the computation of R and the Qk is done, following the 



conclusion of Theorem 3.7 we will compute = {(/%•, Vj) \ j = 1, • • • , n}, derive from the 

Lyapunov exponents lj :— Re(fj,j) the stability of the periodic orbit T and from the eigen 
vectors {vi, . . . , v n } of R we will construct the tangent bundles as defined in Definition 3. 
and given by the formula (|43| 



Computation of R and Qk 
To begin with, let us explicitly write the Jacobian 



and, as consequence, the coefficients Ak 



An 



p-e -i 
e e 





-a 



A k 



-f 3 












-u 1 




-P. 













-a 


a 






jfe > 1. 



The hypothesis (46) for ^ to lie in a ball centered at implies that ||-4.|| s * < oo. Although 
this bound is sufficient to proceed with the computational process, we want to stress out 
that precise informations are known about the |*4fc|oo of the tail elements of the sequence 
Indeed it can be easily seen that 



\Ak\oo < V2r 



1 



VJfc > My 



(48) 



The computation of the approximate solution R, Qk,i, Qk,2 has been addressed as fol- 
low: consider the approximation ^(t) — Y^\k\<M ^fce lfe27rt / Tl of the periodic orbit ^(t) and 
numerically solve system (47) up to time 2f 7 . Denote by y(2f 7 ) the obtained result and 
numerically compute 

TZ = log(y(2f 7 )). 

Neglect the imaginary part and consider only the real part. Then numerically integrate 
the system ^ up to time 2f 7 with 1Z in place of R yielding the solution Q(tj). Fix the 
positive finite dimensional parameter m and compute from Q(tj) the matrices Qk,i, Qk,2, 
respectively the real and imaginary part of the Fourier coefficients with \k\ < m. Finally the 
vector (1Z, Qk) is considered as starting point for a Newton iteration scheme applied on the 
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finite dimensional reduction defined generally in (111. Denote the output of the iterative 
process by x = (R,Qk), that is an approximated solution f^ m \x) w up to a desired 
accuracy, where f^ 1 is defined in (11). 



Consider Af- given by (13). Note that in the case of the three-dimensional vector field 



(44), Afe is a 6 x 6 matrix and one could compute its inverse analytically using the mathe- 
matical software Maple. After having computed A^" 1 one needs to check that the chosen m 
satisfies m > K where K is the same as in Lemma |2.2[ otherwise increase m. 

Then for a choice of M > m and 2 < s < s* one can compute the coefficients Yk 
and Zk, k = 0, . . . , M and Ym-, Zm as shown in Section 2.5 It only remains to define the 
computational parameters Lk introduce in ( |33[ ). In the computation presented here Lk has 
been chosen as 



L k = max{M + M 7 } + k. 



(49) 



This choice assures that the tail elements Hq , Hk in ( 34 ) only contain the terms Aj satisfying 
| ^4. j | oo < y/2r 7 wj s . Therefore the subsequent estimate for hk can be improved by replacing 
||*4|| s * with r 7 . 

Again, the knowledge of the particular behavior of the coefficients Ak allows to provide 
a better estimate for Zm- Indeed note that |i?e(.4fc)| < cw \Re(Ak)\ + w^ s 1„, where Ak 
denotes the matrix Ak with the entries £ in place of £ and the same holds for \Im(Ak)\- 

Therefore \Re(A k )\oo + \Im{A k )\oo < V^|6fe|oo 
|fc| > M 7 . 

Therefore, the computation of the bound for |c£ ijoo when k > M, necessary for the 
definition of Zm, has been slightly modified as follows. 



for 1 < \k\ < M 7 and (481 for 



|Cfc,l| 



J2 (\Re{A 3 )\ + \I m {A 3 )\)w7 s t n 



i+j=k 
\l\jtk 



J2 (1^(^)1 + \Im(Ai)\)wf'l n 



i+j=k 

&0,2k 



< (\ Re ( A i)\ + \Im(A j )\)wl~ s l n + 2r 7 w7 s w7 s %J: 



i+j=k 

|j|<M 7 



l+j = k 



Then, passing to the infinity absolute value, for any k > M 



3=1 



2nr y w j Sw i 
i+j=k 
\l\jtk 




M 



^ I s Af s -!(s- 1) 



M 

1 + *E* 



3 = 1 



I s M s -!(s - 1) 



(50) 



Vm 



(51) 



Computational results 

For the choice a = 10, j3 = 8/3 it is known that there exists a branch of periodic solutions 



parametrized by p joining a Hopf bifurcation at p 
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24.736 and a homoclinic point 



at p Ri 13.9265. Figure [lja-b) shows the bifurcation graph and some of the periodic orbit 
of the continuous family. 
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(a) (b) (c) 

Figure 1: (a) The partial bifurcation diagram for the Lorenz system. The labelled points correspond 
to the values pi. (b) Some of the periodic orbits on the family joining the Hopf Bifurcation and the 
homoclinic point, (c) The periodic solutions corresponding to p — pt. 



The computation of the rigorous enclosure of the periodic orbits and successively of their 
tangent bundles have been performed for a set of different periodic orbits of the Lorenz 
system lying on the mentioned bifurcation branch and corresponding to values of p 

pi = 18.0815, p 2 = 18.6815 



p 3 = 20.8815, pi = 23.8815, p 5 = 24.1816 



Figure 1(c) reports the graphics of the numerical approximation ji of the orbits correspond- 
ing to the choice pi, while the Table ^ contains the computational parameter M 1 that have 
been chosen for the rigorous enclosure of the orbit, the period f 7 and the radius r 7 resulting 
from the computations. The growth rate s* has been fixed s* = 2 for all the cases. 



# sol 


M 7 


T 7 


r 1 


1 


32 


1.027854 


6.844864508150837e 


-09 


2 


30 


0.978271 


7.151582969846857e 


-09 


3 


26 


0.822883 


4.260379031142465e 


-09 


4 


20 


0.683813 


5.368959115576269e 


-09 


5 


30 


0.672595 


2.360935240171144e 


-08 



Table 1: The i-th row concerns the computation of the periodic orbit for the Lorenz system corre- 
sponding to p — pi. M 7 is the finite dimensional reduction parameter chosen in the computation, 
r 7 the approximated period of the solution and r 7 the resulting enclosing radius. 

In the Appendix the first 15 Fourier coefficients of 71 and 74 are listed. As shown in 
Figure[2j one can notice that the Fourier coefficients of the five orbits under consideration are 
decaying to zero with a different speed. This is due to the fact that the closer we are to the 
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1 23456789 10 

k 



Figure 2: Norm of the Fourier coefficients of each of the periodic solutions 7i 

homoclinic orbit, the flatter the periodic solution is, meaning that a larger number of Fourier 
coefficients contributes to the Fourier expansion, hence leading to a slower decay. Table [2] 
contains information about the computation of the fundamental matrix solution associated 
to each of the previous periodic orbits. More precisely, it contains the finite dimensional 
reduction parameter m, the computational parameter M and the resulting radius r. 



# sol 


m 


M 


r 


1 


100 


180 


1.98645943e- 


05 


2 


90 


140 


7.52145121e- 


06 


3 


80 


80 


9.66152623e - 


07 


4 


60 


66 


9.91268997e - 


07 


5 


60 


70 


3.77687574e - 


06 



Table 2: Computing the fundamental matrix solution for each of the periodic orbit 7 4 . m and M 
are respectively the finite dimensional reduction parameter and the computational parameter that 
have been chosen, r is the radius of the ball centered approximate solution in Q" within which a 
genuine solution of Q exists. 

Some of the radii polynomials Pk(r) built during the computation of solution #4 have 
been plotted in Figure [3] The bold line on the x-axis remarks the interval 

INT = [9.91268997 • 10" 7 , 1.4574858482 ■ 10" 3 ] 

where all the radii polynomials are negative. 

From the computations we noticed that the odd Fourier coefficients of Q(t) are almost 
vanishing, suggesting that Q(t) is a r 7 periodic function, rather than 2r 7 periodic. This is 
not in contradiction with Floquet Theorem. Again in the Appendix we report the numerical 
approximation R and the first even Fourier coefficients Qk for the solution #1 and solution 
#4. As in the previous case, the Fourier coefficients Qk corresponding to periodic orbits 
closer to homoclinic decrease slower. This justifies the fact that larger values of m and M 
were necessary to obtain successful computations. 
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0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 " s lu 15 

' X10" 4 

Figure 3: Plot of some of the radii polynomials Pk(r) constructed in the computation of the 
fundamental matrix solution associated to 74. On the right: magnification close to r = 0. The red 
line denotes the interval INT where all the Pk (r) are negative. 



We now have all the ingredients necessary to construct the tangent bundles: first we 
compute the intervals containing the spectrum and the eigenvectors of each the interval 
value matrix R, then, in light of Theorem |3.7[ the multiplication of the stable and unstable 
directions with the function Q(9) yields the tube enclosing the complete stable and unstable 
bundles. 

Table[3]lists the Lyapunov exponents of the periodic orbits, as defined in Definition 3.6 



and it also contains the radius of the intervals enclosing the stable and unstable eigen-couple 
of R while in Figure [4] the tangent bundles are depicted. In Appendix the complete list of 
the eigen-decomposition of the interval matrices R is also provided. 



Sol# 


Center Radius 


1 


-14.2953855130260 6.801248614 • 10" 5 
0.6287188463595 9.510853040 • 10" 4 


2 


-14.2174898849454 2.814282939 • 10" 5 
0.5508232182790 3.733964559 ■ 10~ 4 


3 


-13.9620493680589 2.774785811 • 10" e 
0.2953827013923 3.653168667 • lO" 5 


4 


-13.7210150091049 2.544262339 • 10~ 6 
0.0543483424385 5.248456341 • 10~ 5 


5 


-13.7013292393391 9.720262854 ■ 10" e 
0.0346625726730 3.336819199 • 10~ 4 



Table 3: Lyapunov exponents for each of the periodic orbit ji. For each solution we report the 
center and the radius of the interval vectors enclosing the exponents. Note that we could prove the 
existence of the eigenvectors Vj associated to within the same accuracy given by r. 



4.2 £ 3 -model: non orientable tangent bundles 

It is known that if a Floquet multipliers of a periodic orbit is negative, then the corresponding 
tangent bundle is not orientable. Moreover, in the case of a saddle periodic orbit of a three- 
dimensional system, the two non-trivial Floquet multipliers are real and their product is 
positive. Therefore both the tangent bundles are either orientable or not orientable and, in 
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Figure 4: Plot of the tangent stable (turquoise) and unstable (red) bundles of each of the 
periodic orbits ji. The central figure concerns Sol# 2, with a magnification, while figures 
(a)-(b)-(d)-(e) concern respectively 71, 73,74,75- 



2G 



the latter case, they are topologically equivalent to a Mobius strip, see [55] . 

An example of a dynamical system with periodic orbits that exhibit this behavior is the 
so called £ 3 -model considered in 



x 

V - 



u 

z 

ax 



(3y- 



(52) 



For fj = 2, as a varies, the periodic orbits of system (52 1 produce an interesting bifurcation 
diagram. We refer to [35] and |37j for a detailed analysis of the bifurcation diagram and 
on the genesis of periodic orbits, called twisted periodic orbit, with non orientable invariant 
manifolds. We focus on a particular twisted periodic orbit corresponding to a — 3.372 lying 
on the branch emanating from a period-doubling bifurcation that occurs at a w 3.125. 

Following the same procedure as before, we rigorously compute the enclosure of the 
periodic orbit "f(t) and subsequently the enclosure of the matrix R and of the matrix function 
Q(t), hence producing an explicit Floquet normal form as in <[2). Then, we extract the 
necessary stability parameters and we recover the stable and unstable tangent bundles using 
(43). Figure [5] shows the resulting bundles. 

Having computed the intervals enclosing the period r of the orbit and the eigenvalues of 
R, we realize that the absolute values of the two nontrivial Floquet multipliers satisfy 



I cri I <= [7.037235782193 ■ 1CT 3 
I £7 2 1 e 1.526609276443494 



7.037944324307 • 10~ 3 ] (= A st ) 
1.528421395487018] (= A„ st ) 



To conclude we emphasize the role played by the continuous function Q{9) in the con- 
struction of the tangent bundles. As proved in Theorem |3.7[ as 9 changes, the eigenvec- 



Q(0)vj, 



where 



tor Wj of $e(T) associated to the Floquet multiplier £7j is given by w 
Vj is the eigenvector of R relative to the eigenvalue fij. The function Q{9) is continu- 
ous and 2-r-periodic, but the tangent bundles are smooth manifolds. That implies that 
wj = Q(r)vj has to be an eigenvector of $(r) associated to the Floquet multiplier (jj, i.e. 
span{vj} = span{u>j}. In the case of the Lorenz system Q(t) turns to be the identity 
matrix, therefore the last relation is simply verified. But in case of the £ 3 -model and in 
general when the bundles are not orientable Q{t) need not be the identity matrix. Indeed, 
in the considered example, Q{t) results to stay in a small interval around 



Q 



-1.663148705259924 
1.323227706784005 
0.936264065136750 



-1.018593776943882 
1.032472501143805 
1.438097884773345 



-0.446703151142258 
0.891338521225762 
-0.369323795883880 



Denoting by R, r, y\ the centers of the intervals the genuine R, r and v\ belong to and 
defining $ = Q(r)e RT the numerical approximation of $(T), we compute 



<5>Qv 



-0.002642211417990 
-0.003294408641461 
0.011434535907588 



Qv 



0.375442644619642 
0.468116019931565 
-1.624780050494352 



The component-wise ratio between the two computed vectors is a% = —7.037590044326 • 
10~ 3 ±10~ 13 , whose absolute value is indeed in the interior of A st . If the unstable eigenvector 
t> 2 is considered, the same operations produce cr 2 = — 1.527515067244305 ± 10~ 13 . Although 
not rigorous, these computations confirm the above theoretical discussion and moreover 
provide a method to recover the sign of the Floquet multipliers, information that is not 
possible to achieve following the presented computational technique. 



27 



(a) 




5 Acknowledgments 



We would like to thank Marcio Gameiro and Jason D. Mireles James for helpful discussions. 



6 Appendix 



Period and Fourier coefficients of 71 and 74. 



Solution # 1 

T 7 = 1.027854840752128, f = 



-4.606354666884038 
-4.606354666884038 
13.533127936581090 



(3 



-2.457589444025310 - 0.734232230617879s 
-2.008759805775985 - 2.236534819812307s 
3.644641521462053 - 2.592383847330301s 

«4 

-0.066251105459624 - .023092599715315s 
-0.009785901431526 - . 185087447977630; 
0.181193087615871 - 0.007305965402112s 



-0.001215481367863 - 0.000664841070850i 
0.001629398704351 - .005865931787398i 
0.005830102115005 + 0.001711097840891; 

iio 

-0.000021217174006 - 0.000017849897868i 
0.000087897665531 - 0.000147548597224i 
0.000146654700161 + 0.000089133299138i 



(13 



-0.000000344798445 - .000000437771341s 
0.000003134076053 - 0.000003177810272s 
0.000003154570627 + 0.000003153649381i 



-0.840331595816763 - 0.280556868231764s 
-0.497327754727251 - 1.307931343042055s 
1.373048803380924 - 0.590924707028489s 



-0.017586386797219 - 0.007027398326009s 
0.003892543961801 - 0.060779408215756; 
0.059967166034911 + 0.005022757182874s 

«8 

-0.000317263424301 - 0.000201467445482s 
0.000667978649296 - 0.001752989942270i 
0.001743196289261 + 0.000687999808866i 

111 

-0.000005425393706 - 0.000005230883352s 
0.000029748123926 - 0.000041712327941s 
0.000041443384209 + 0.000030059169428; 

|l4 

-0.000000085406182 - 0.000000125225422s 
0.000000986282445 - 0.000000856137980s 
0.000000849430424 + 0.000000991 127939i 



-0.244051806246999 - 0.079599377946577s 
-0.098076628971404 - 0.527159479549944s 
0.516674965708581 - 0.103578564294281s 



_S6_ 



-0.004633846567037 - 0.002168050944621s 
0.003318015096906 - 0.019163826327315s 
0.019005700003600 + 0.003642138482031s 

«9 

-0.000082326458994 - 0.000060305586234s 
0.000249451961957 - 0.000513234481766s 
0.000510298663577 + 0.000254394262682s 

£l2 

-0.000001374889998 - 0.000001519316512s 
0.000009770046204 - 0.000011604812976s 
0.000011525197547 + 0.000009848310635s 



iis 



-0.000000020840071 - 0.000000035583174s 
0.000000305435097 - 0.000000226673423s 
0.000000224755146 + 0.000000306618231s 



Solution # 4 

T 7 = 0.683813590045746, f 



-7.521252250993276 
-7.521252250993276 
22.399077327399255 



ii 



S2 



_la_ 



-1.246453092091490 - 1.262394967959499s 
-0.117098081743200 4- 1.579694736267260s 
-1.138858133842916 - 1.528289275687748s 



-0.114587819240451 - 0.194828460289262s 
0.015748373138453 + 0.188342422139680s 
-0.143528439361466 - 0.149121347835789s 



-0.008589210614442 - 0.020708730881492s 
0.002967722490543 + 0. 020459531930793s 
-0.016769839350022 - 0.016159568619679s 



(5 



-0.000619073947482 - 0.001950094567689i 
0-000325236831105 + 0.001967969477674s 
-0.001814442573257 - 0.001764257857272s 



-0.000044052745410 - 0.000171023517657s 
0.000031364652698 + 0.000173341352082s 
-0.000188148978551 - 0.000185385380347s 



-0.000003073584415 - 0.000014085485779s 
0.000002859395743 + 0. 000014272552631s 
-0.000018837629098 - 0.000018676999462s 



£7 



S8 



(9 



-0.000000208407126 - 0.000001089211660s 
0.000000251244910 + 0.000001104010427s 
-0.000001824392597 - .000001813988748i 



-0.000000013587818 - 0.000000078416175s 
0.000000021464593 + 0.000000079624828s 
-0.000000171368842 - 0.000000170686207s 



-0.000000000836854 - 0.000000005128148s 
0.000000001792309 + 0.000000005227527s 
-0.000000015658534 - 0.000000015616179s 



ilO 

-0.000000000046957 - 0. 000000000284715s 
0-000000000146750 + 0.000000000292804s 
-0.000000001395354 - 0.000000001392990s 



Cll 

-0.000000000002182 - 0.000000000010252s 
0.000000000011806 + 0.000000000010900s 
-0.000000000121501 - 0.000000000121396s 



iu 

-0.000000000000052 + 0.000000000000357s 
0.000000000000935 - 0.000000000000305s 
-0.000000000010352 - 0.000000000010351s 



0.000000000000006 + 0.0000 00000135s 
0.000000000000073 - 0.000000000000131s 
-0.000000000000863 - 0.000000000000865s 



0.000000000000002 + 0.000000000000021s 
0.000000000000006 - 0.000000000000020s 
-0.000000000000070 - 0.000000000000071s 



0.000000000000000 + 0.000000000000003s 
0.000000000000000 - 0.000000000000003s 
-0-000000000000005 - 0.000000000000006s 



Numerical approximation R and even Fourier coefficients Qk- 



Solution # 1 



-10.508958375451483 
1.367770562612481 
-6.918853545877750 



6.244108010218356 
5.059391467543374 
5.863201994753524 



-7.445538972862637 
-10.140640221489871 
-8.217099758758689 



Qo ■ 



1.411735844583484 
-1.110555333911319 
0.843423174876108 



-0.999238898309471 
0.141112697583195 
0.200767676367087 



1.303728854375973 
0.194620640182095 
-0.730409462001914 



Q2 
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-0.263681129667594 - 0.434148218933761; 

0.495149611356963 + 0.020091337974708; 
-0.125150133979381 + 0.155543962133145; 



0.245502648963791 + 0.051118262858919; 
0.018943994759567 - 0- 173000962596422i 
0.072778675242146 4- 0.276759663116221; 



-0.264146740709986 - 0.219962925571520; 
0.034736062134259 + 0.070721806450736; 
0.234684069057806 - 0.137872983068902; 



0.003223434126082 -f 0.072191532567703; 
0.088047497733925 + 0.044110570228044; 
-0.138910366811394 4- 0.089747331219474; 



0.032229467843581 4- 0.043355074683078; 
-0.000475098772244 4- 0.078582831969487; 
-0.085890896774482 4- 0.011537815519094; 



0.014837677547204 + 0.016791777242045; 
-0.011244165701100 + 0.044396287876498; 
-0.042670724306854 - 0.009276118631273; 



0.005198646646053 + .006059114013686; 
-0.008342719368296 + 0.019119854569593; 
-0.018320486491551 - 0.008387560951806; 



0.001642111476789 + 0.002113215989109; 
-0.004468250361663 + 0.007266315842478; 
-0.007054400616465 - 0.004581594079138; 



0.000489747959995 + 0.000714963537264; 
-0.002038899547182 + 0.002555634111674; 
-0.002503668907675 - 0.002080001556342; 

Solution # 4 



R = 



Qo = 



0.043888342026888 - 0.018225551633217; 
0.184261504574051 + 0. 121001510966144; 
-0.196713672633492 4- 0.151158740280049; 



0.020354746043569 + 0.006072390169069; 
0.019156656540045 + 0.040757474867593; 
-0.042700076948069 + 0.023165829211282; 



0.002761640132931 4- 0.000187559216871; 
0.003017836182538 + 0.007295810537329; 
-0.007076735073531 + 0.002990041055165; 



Q4 = 

0.154600837185400 - 0.083425433638396; 
0.205199858834844 4- 0.061631720475212; 
-0.038756214412790 + 0.227619738561727; 

Qe = 

0.063509812104330 - 0.036923452979636; 
0.122476382683268 + 0.065280023891 177i 
-0.060779388764130 + 0- 1 18929071081753i 

Qs = 

0.023655160613775 - .012807406159434; 
0.053040611505731 + 0.039952372292969i 
-0.039817529715657 + 0.050937744964250i 

OlO = 

0.008265722985980 - .003957671 147668» 
0.019900084741206 + 0.019629584406457i 
-0.019826346511011 + 0.019224271692316i 

<5l2 = 

0.002770838660308 - .001 130105417024; 
0.006811683852990 + 0.008504503427600i 
-0.008590105865316 + 0.006619388137483; 

Ql4 = 

0.000900858813617 - 0.000302058076912; 
0.002173495066311 + 0.003390709655169; 
-0.003415963753540 + 0.002118659036427; 



Q2 = 

0.237009816130334 - 0. 161848354458699i 
0.327547824090844 - 0.006620885935044; 
0.033125259013010 + 0.469308310958510; 

«4 = 

0.031468524011727 - 0.044530642983650; 
0.108300555430671 + 0.012024500455431; 
-0.006988160820195 + 0.106913456324494; 

<56 = 

0.002550968100164 - 0.005879343921945; 
0.018305847192402 + 0.001100730065600; 
-0.001327830255369 + 0-017843392559588; 



-0.224496433432301 - 0.070860905176751; 
-0.101564118363285 - 0.257855508650553; 
0.277112301560736 - 0.128593935735273; 



-0.104043649831591 - 0.038640713988921; 
-0.036087360927645 - 0.192529574267658; 
0.189091902999670 - 0.038256089963182; 



-0.039506811034992 - 0.015496956153149; 
-0.004186361886575 - 0.098878456196623; 
0.096680260325811 - 0.002950867660488; 



-0.013457029500778 - 0.005897110459635; 
0.003372357081277 - 0.042680213775582; 
0.042021880260296 + 0.004167728643705; 



-0.004317304940509 - 0.002176073388202; 
0.003179105256374 - 0.016650157014229; 
0.016483705220024 + 0.003472785816957; 



-0.001332830306602 - 0.000777133607668; 
0.001809581469994 - 0.006067809090768; 
0.006023316397811 + 0.001899418469947; 



-0.187509342015513 - 0.182145921801054; 
0.028142805119932 - 0.301161573224061; 
0.383336675401313 - 0.015987641512644; 



-0.037967521528517 - 0.020350032723162; 
0.003607621081908 - 0.088848037842172; 
0.087287750910915 + 0.001272272562888; 



—0.004808213884619 - 0.001470031764531; 
-0.000373601315687 - 0.014656539554954; 
0.014369458507765 - 0.000147876219115; 



-10.103827000749006 
2.108771563239242 
-6.292527840128125 

0.865350358013670 
-0.413891831683466 
0.495177534990071 



5.011512150268070 
-0.623931925418962 
3.486887629270139 

-0.542407880588934 
0.086095506914414 
-0.049624165086980 



-4.181592133228406 ] 

0.486619976897008 
-2.938907740498710 J 

0.461699549706412 
-0.062238564323011 
0.025622831143080 



0.000289533560511 - 0.000023944541333; 
0.000441108794910 + 0.000976089854840; 
-0.000969488479052 + 0.000417372380704; 



0.027798782570345 - 0.005964538385053; 
0.060031396019271 + 0.114573660552919; 
-0.115071668553051 + 0.057995199816582; 



0.025222741873242 - 0.008833955070603; 
0.077200510619449 + 0.122891998519810; 
-0.123626123118424 4- 0.075719660830317; 



0.021844676611758 - 0.010885662915616; 
0.093733134787201 + 0.122596846519192; 
-0.123324223393918 + 0.092568245260174; 



0.000167141543813 - 0.000637470930420; 
0.002466959960129 - 0.000020947005419; 
-0.000027259876302 + 0.002448108770303; 

Q 10 = l ie - 03* 

0.007649736579992 - 0.063378060920900; 
0.294930056067831 - 0.027358497808776; 
0.023085665062373 + 0.294833391391544; 

<5l2 = Lie - 04* 
-0.000684171241317 - 0.059552377024599; 
0.324298385293782 - 0.061987569005484; 
0.058531951875591 + 0.324753469250377; 

Q 14 = l.le - 05* 

-0.007613608020911 - 0.053546469056524; 
0.334057192499540 - 0.100903036175184; 
0.097987409026114 + 0.334717932778523; 



-0.000509911668000 - 0.000073299205025; 
-0.000203588657299 - 0.001946883870274; 
0.001937426651095 - 0.000166408614554; 



-0.049748597822010 - 0.000258862545890; 
-0.045190737024386 - 0.229328952532320; 
0.229589287166289 - 0.041943841171559; 



-0.045934587657163 + 0.005873935636186; 
-0.075398660882460 - 0.248296827479228; 
0.248928813781281 - 0.072775716802859; 



-0.040603455774486 + 0.010691275157818; 
-0.106940875623844 - 0.251662436543305; 
0.252425413480351 - 0.104748352612175; 



Enclosure of the spectrum and eigenvectors of R. 



Solution # 1 



E .valut 
E .vecto 



Stable | 



Unstable 



-14.295385513026014 
-1.304013849063401 
-0.455394137737842 
-1.045066534133045 



-0.000000000000342 
0.330244752093107 
1.503215970221508 
0.794531403146519 



0.628718846359581 
0.376869068070140 
1.529903371837170 
0.719281153912157 



0.068012486147407 | 0.933273952985148 | 0.951085304085387 
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Solution # 2 



E. values 
E. vectors 



_L 



_L 



-14.217489884945432 | -0.000000000000372 | 0.550823218279069 



-1.311412833044274 
-0.418001967490342 
-1.051413684760184 

0.028142829394895 



0.313018700843905 
1.491825566190409 
0.822481472729119 

0.366070961579701 



0.362984001981858 
1.522008146749871 
0.742787867114333 

0.3733964559921 



Stable 



E. values 
E .vectors 
Had 10~ 4 - 



-13-962049368058929 | 
1.347327907101522 I 
0.192884294913700 
1.071215739018561 



-0.000000000000126 
0.210153254038267 
1.398977370916435 
0.999348750677595 



0.295382701392358 

-0.271285496065970 
-1.447969583732432 
-0.910927145390874 



0.027747858117877 | 0.357120423613300 | 0.365316866760002 



Solution # 5 



E.valw 
E.vecto 
Had 10~' 



E . values 
E .vectors 



Stable | 



Unstable 



-13.721015009104903 | 

1.439298428490600 
-0.266150110668696 
0.926058395748094 



-0.000000000000309 
0.051023279023768 
1.168676456309931 
1.277337843119246 



0.054348342438550 
-0.128543563503133 
-1.252767532277167 
-1.189138369725781 



0.025442623394767 | 0.486579977382052 | 0-524845634189975 



Stabh 



_1_ 



Unstable 



-13.701329239339196 | 
1.451175793211715 

-0-333782791906082 
0.884690830190827 



-0.000000000000487 
0.020133220407690 
1.115780576244332 
1.324623855708436 



0.034662572673008 
-0.100604795866982 
-1.206723375086350 
-1.238425359506485 



0.009720262854618 | 0.336965856255784 | 0.333681919907764 
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